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ABSTRACT 

Context. Young terrestrial planets, when they are still embedded in a circumstellar disk, accumulate an atmosphere of nebula gas. The 
evolution and eventual evaporation of the protoplanetary disk affect the structure and dynamics of the planetary atmosphere. These 
processes, combined with other mass loss mechanisms, such as thermal escape driven by extreme ultraviolet and soft X-ray radiation 
(XUV) from the young host star, determine how much of the primary atmosphere, if anything at all, survives into later stages of 
planetary evolution. 

Aims. Our aim is to explore the structure and the dynamic outflow processes of nebula-accreted atmospheres in dependency on 
changes in the planetary environment. 

Methods. We integrate stationary hydrostatic models and perform time-dependent dynamical simulations to investigate the effect of a 
changing nebula environment on the atmospheric structure and the timescales on which the protoatmosphere reacts to these changes. 
Results. We find that the behavior of the atmospheres strongly depends on the mass of the planetary core. For planets of about Mars- 
mass the atmospheric structure, and in particular the atmospheric mass, changes drastically and on very short timescales whereas 
atmospheres around higher mass planets are much more robust and inert. 
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1. Introduction 


At early stages in the evolution of planetary systems, newly 
Formed protoplanetary cores composed of solid material orbit the 
host star still embedded in the circumstellar disk. These cores in¬ 
teract gravitationally with the disk gas and thus accumulate an at¬ 
mosphere, which at its outer margins blends into the surrounding 
disk nebula environment. Eor more massive cores, more substan¬ 
tial atmospheres form accordingly, and, if one follows the core¬ 
accre t ion paradigm of gas-planet formation dPerri & CameronI 
1971 Im izunollT9^ . this finally leads, for sufficiently massive 
:ores, to runaway accretion and the subsequent formation of a 
gas giant. 

Eor the numerical study of embedded planetary atmospheres, 
it is necessary to adopt some assumptions and simplihcations. 
These include; 


1) Spherical symmetry: Eirst, the Hill radius defined as 


^Hill - 


3mJ 


a , 


( 1 ) 


with a the semimajor axis of the planet, gives an approximation 
for the sphere where the gravitational potential is dominated by 
the planet with a mass Mpi over the host star with a mass M*. 
Another limit comes from the requirement that the disk medium 
is essentially uniform along the perimeter of the planetary at¬ 
mosphere. The vertical structure of the protoplanetary disk can 
be approximated by a hydrostatic stratification with a pressure 
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scale height Hp. Assuming a vertically isothermal disk with a 
temperature T, the pressure scale height can be expressed as 


RTa^ 
]j GM 


( 2 ) 


where K is the specific gas constant, and G is the gravitational 
constant. Combining these two arguments, spherical symmetry 
is a good approximation as long as r and r «; T/p. In this 

paper, we always assume that £ Hp. Spherical symmetry 
also implies that rotational momentum of the atmosphere with 
respect to the planetary core can be neglected. 


2) Stationary equilibrium: Dynamical, i.e. hydrostatic, equi¬ 
librium, is certainly well fulfilled during most phases of atmo¬ 
spheric evolution. Thermal equilibrium, however, is much harder 
to justify and only adopted pragmatically to make the problem 
numerically tractable. The assumption of thermal equilibrium 
implies that the luminosity in the atmosphere is radially constant, 
therefore the luminosity at the surface of the core is established 
as a m ain constituting parameter. Traditionally, (iHavashi et al.l 
Il979l) this luminosity has been associated with the energy re¬ 
leased by accreted planetesimals (assumed as continuous stream) 
as they impact on the surface of the core. Equating the released 
energy with the gravitational energy of the stream of planetesi¬ 
mals, the planetary luminosity Lp\ is directly related to the accre¬ 
tion rate Macc of the planetary core 


Tpi — GAlp\]M^QQ 




(3) 
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where Mpi and Rpi are the mass and the radius of the plane¬ 
tary core, respectively. The planetesimals are assumed to travel 
through the atmosphere without energy losses so that all of the 
gained gravitational energy is released on the planetary surface. 

Considering these two assumptions, the structure of a plane¬ 
tary atmosphere is described by the usual stellar structure equa¬ 
tions that can be solved readily by established numerical meth¬ 
ods. Such hydrostatic and spherical symmetric models of pri¬ 
mary hydroge n atmospheres around terr estrial planets have been 
calculated bv iHavashi et al.l (Il979ll and iNakazawa et ^ (Il985h 
using a constant and an an alytically prescribed du st opacity, 
respectively. More recently, llkoma & Gendal (l2006h improved 
upon these earlier results with the adoption of realistic equa¬ 
tion of state and opacities. Particularly important for the atmo¬ 
spheric structure is detailed dust opacity data, as the evapora¬ 
tion/condensation of various dust species directly manifests in a 
sequence of atmospheric dust layers. 

This work deals with two distinct numerical methods: sta¬ 
tionary modeling and time-dependent simulations. After intro¬ 
ducing the physical scenario in Sect. o the method and results 
of stationary modeling are presented in Sect. |2] Then Sect. [3 
describes the time-dependent scheme and the results of the dy¬ 
namic simulations. Finally, in Sect.|4] we discuss both the effect 
of numerical parameters and the astrophysical relevance of the 
results as well as the implications for the evolution of habitable 
planets. 


1.1. Physical scenario and boundary conditions 


From statistics of IR-excess observations (iHaisch et al.l l200il: 
iHillenbrandl l2005h it is readily evident that the circumstellar 
disks have a typical lifetime of a few Myr. After about 10 Myr 
almost no protoplanetary disks are left. The mechanism of disk 
evaporation is still not well understoo d, though some plausi¬ 
ble scenarios have been devised (e.g.. iHollenbach et'aL 199^ 


iJohnstone et al.ll9^IClarke et al.l2001blScallv & Clarke 200 Ih . 

It is, however, clear that the disappearance of the nebula environ¬ 
ment will have significant consequences for the formerly embed¬ 
ded planetary atmospheres. 

In this study, we focus on the direct hydrodynamical effect 
of the disk evaporation on the atmosphere, ignoring many poten¬ 
tially important physical processes resulting from the exposure 
of the planetary atmosphere to the radiation field of the host star, 
such as escape from a thermosphere heated by the high extreme 
ultraviolet and soft X-ray luminosity (XUV) of young stars 


(lYelle 


l2013FlLammer et al.il2013i 12014 . The interaction of an evapo 
rating disk with the planetary atmosphere around a cooling plan - 
etary core recently has been studied by llkoma & Horil ( 20121) . 
Whereas this study considers the dynamic reaction of atmo¬ 
spheres, which initially are bo th in dynamic and thermal equi¬ 
librium, llkoma & Horil (1201 2l) use quasi-static calculations and 
cover the atmospheric evolution as the atmosphere and the solid 
core cool down from an initially hot state. 

In terms of hydrostatic balance, the disk nebula environment 
establishes a fixed gas pressure at the top of the planetary atmo¬ 
sphere. Once the circumstellar disk vanishes, this gas pressure 
cap is released to a large degree. In our models in particular, we 
reduce the gas density at the outer boundary by 15 orders of mag¬ 
nitude. Starting from an initial density of p = 5 x 10“'° g/cm°, 
in accordance w ith estimates of the initial mass planetary nebula 
(iHavashil 19811) at an Earth orbit, the outer boundary density thus 
drops to a value of p = 5 x 10“^^ g/cm°. 


Koskinen et al. 12013^ lLammedl2013t lErkaev et alJ 


Technically, these very low densities are in a range where the 
fluid approximation of hydrodynamics becomes invalid. In our 
calculations, we neglected this effect to represent, in the contin¬ 
uum description, the scenario of p ^ 0 (and thus P ^ 0) at 
the outer margin of the planetary atmosphere. Moreover, in real¬ 
ity the very thin outer region of the atmosphere is dominated by 
stellar wind and stellar irradiation. Consequently, the gas is ion¬ 
ized to a significant degree and a plasma can exhibit collective 
phenomena arising out of mutual interactions of many charged 
particles. This is usually described by the so-called dimension¬ 
less plasma parameter A, which basically characterizes the num¬ 
ber of particles within a Debye sphere, i.e., A - n 47:A^, where 
n is the number density of electrons and To is the Debye length. 
The larger this plasma parameter is (typically A ^ 10° in the 
interplanetary medium), the more collective interactions are tak¬ 
ing place because of the decreasing effective shielding volume. 
Hence, a situation is established that still allows for a hydrody¬ 
namical description of the fluid (e.g.. lChoudhuril 19981) . 

In terms of astrophysical evolution of planetary systems, the 
evaporation of disks is a quite fast process. Statistics of transi¬ 
tional disks indicate a c haracteristic transition timescale of about 
(0.1 - 1) X 10° years (iMuzerolle et al.ll2010l) . As we are inter¬ 
ested in the dynamic reaction of the atmosphere, we have to im¬ 
plement this density change on timescales smaller than the dy¬ 
namical timescale of the atmosphere, independent of the realistic 
evaporation timescale of actual disks. In Sect. 13.21 we show that 
for our atmospheric models the time to adjust to changes in the 
outer boundary density ranges between ~ 10^ and ~ 10*° years. 
Eor consistency and comparability, the evaporation of the disk 
is included in all dynamical simulations on a common timescale 
shorter than this range, in particular, over a period in time of 10° 
s (^ 32 yr). The temperature at the boundary condition is kept 
fixed at 200 K, which is about the radiation equilibrium temper¬ 
ature at an Earth-like orbit at 1 AU around a young solar-type 
star throughout the whole simulation run. 


2. Stationary models 

2.1. Method 


The stationary models are computed from the hydrostatic struc¬ 
ture equations 


dm 


4nr^p , 


dP Gm 
dr ^ ’ 

and 

dT _ dP 
~ ' P 


(4) 

(5) 

(6) 


The meaning of the symbols we used is complied in Tab.[T]These 
equations are solved using the initial model integrator of the 
TAPIR-Code (The adaptive, implicit RHD-Code). The TAPIR- 
Code provides a framework for the implicit solution of the equa¬ 
tions of radiation hydrodynamics on ada ptive grids and pr evi- 
ously has been ap plied to spherical flo ws dPorfi et al.ll2006h . 2- 
dimensio nal flows (IStokl & Dorfll20(y^ . convectio n modeling in 
Cepheids ( St6kj2008l). and planetary atmospheres (lErkaev et al.l 
l20l4lLammer et al.ll20l4 . 

The logarithmic temperature gradient V is either determined 
from radiative transport 


3kkLpiP 

64ncrGmT'^ 


(7) 


Article number, page 2 ofM 



























































A. Stokl et al.: Hydrodynamic simulations of protoatmospheres 


Table 1. Meaning of Symbols. 


Symbol 

Meaning 

r 

t 

Radius from the center of the planet 

Time 

m 

Integrated mass (core mass -f atmospheric mass) 

P 

Gas pressure 

T 

Gas temperature 

P 

Gas density 

e 

Specific gas internal energy 

u 

Gas velocity 

It' 

Convective velocity perturbation 

J 

Zeroth moment of radiation intensity (oc radiation energy) 

H 

First moment of radiation intensity (oc radiation flux) 

K 

Second moment of radiation intensity (oc radiation pressure) 

O 

Gravitational potential 

CT 

Stefan-Boltzmann constant 

G 

Gravitational constant 

K 

Specific gas constant 

cp 

Specific heat at constant pressure 

r 

Adiabatic index 

Cs 

Adiabatic sound speed 

kr 

Total (gas & dust) opacity, Rosseland-mean 

Kp 

Total (gas &. dust) opacity, Planck-mean 

^ d\nT 

^ ~ a\nP 

Logarithmic temperature gradient 

^1 

Planetary luminosity 

Mpi 

Mass of the planetary core 


Mass of the host star 


or, in case of convective instability, calculated from the turbulent 
convection model yielding a temperature gradient Vconv, usually 
close to Vad- In case of the integration of hydrostatic models, we 
use the stationary limit of the time-dependent formulation of the 
turbulent convective energy equation given below in Eq. |2T] 

As noted in the introduction, spherical symmetry is a good 
assumption for r and breaks down completely at 

r > therefore, the natural point to implement the bound- 
ary conditions is the Hill radius. Some previous studies (e.g., 
Ilkoma & Gendal[2006l: iRafikovIbOObI: llkoma & Horill2Ql^ used 

the minimum of Hill radius and Bondi radius (where the lat¬ 
ter typically is about one order of magnitude smaller) for the 
placement of the outer boundary. The Bondi radius relates the 
the gravitational potential to the enthalpy of the gas and thus 
provides an indication to which extent the atmosphere is gravi¬ 
tationally bound to the planet 

GMpi 

^Bondi - - 2 — ’ 

ct 

where Cs is the speed of sound in the nebula. The Bondi radius 
provides a good estimate for the point where the atmospheric 
density structure attaches to the surrounding nebula gas. There 
is, however, no reason not to extend hydrostatic models beyond 
the Bondi radius into the nebula environment u p to the ulti¬ 
mate limit of spherical modeling at the Hill radius (IMizuno et al.l 
Il978h . In view of the dynamic calculations in Sect. [2 where it 
is instrumental to include a domain beyond the immediate plan¬ 
etary envelope structure, we placed the outer boundary of our 
models on the Hill radius throughout. For the calculation of the 
Hill radius, we adopted a — lAU and M* = IMq as standard 
parameters. 

The accretion rate of planetesimals, and in consequence the 
planetary luminosity, is unfortunately a highly uncertain param¬ 
eter. Almost the only indication comes from the observed disk 
lifetime of some 10® years and the fact that the planetary cores 
must have formed within that period of time. Accordingly, we 
adopted a mass accretion rate Macc/ATpi of lO^^yr^' as the stan¬ 
dard value. 

In addition to accretion, the models also allow for radio¬ 
genic heat production, which is included proportional to the core 



Fig. 1. Density as a function of radius in protoplanetary atmospheres 
around cores of 0.1 and 1 and for nebula gas densities varying 
from 5 X 10“'® g/cm^ to 5 x 10“^^ g/cm^ with steps of one magnitude. 
The inner and outer boundaries are defined by the radius of the solid 
core and by the Hill radius, respectively. 


mass using a reference value of 10^* erg s“* Mg,* for the young 
Earth (IStacev & Davisll2008h . This contribution, however, never 
becomes relevant as the accretion luminosity is more than two 
orders of magnitude larger for all models calculated. 

For the atmospheric gas, we assumed a homogeneous chem¬ 
ical composi tion with solar abundan ces. The data for gas opacity 
is taken from iFreedman et ak| ( 2()()8|) and for th e equation of state 
we adopted the results of Saumon et al.l (Il995h . Dust is included 
in the model with a fixed dust grain depletion fa ctor / using dust 
opacity coefficients from ISemenov et ^ (l2003h . The dust deple¬ 
tion factor / specifies the amount of dust in the atmosphere rel¬ 
ative to the local dust condensation/evaporation equilibrium. As 
dust opacities dominate over gas opacities by many magnitudes 
in large parts of the atmosphere, / in effect scales the overall 
opaqueness of the atmosphere. To allow for dust depletion pro¬ 
cesses such as rain-out in the atmosphere and dust segregation 
within the circumstellar disk structure, / is usually set to values 
/ < 1 and we adopted a standard value of / = 0.01. 


2.2. Results 

Figure[T]illustrates the dependence of the hydrostatic atmosphere 
structures on density variations at the outer boundary for cores of 
0.1 M® and 1 M®. The very different behavior in the two model 
sequences is striking: whereas the atmosphere around the Mars- 
mass core almost completely escapes with the removal of the 
surrounding nebula gas, the structure of the atmosphere around 
the Earth-mass only changes gradually and substantial parts of 
the atmo sphere remain gravitatio nally bound. This confirms the 
results of llkoma & Gendal(l2006h . who found an almost linear re¬ 
lation between atmospheric mass and nebula density for a Mars- 
mass core in contrast to a much less dependent atmospheric 
mass for an Earth-mass planet. In terms of a bsolute numbers, 
our re sults are in agreement with those from llkoma & Gendal 
(1200^. The ma i n phy sical difference between this study and 
llkoma & Gendal (l2006h is the data for dust and gas opacities. 
Other numerical differences, such as the treatment of convective 
transport and the placement of the outer boundary condition, are 
less important. 

Figure |2] illustrates the disparity in atmospheric behavior for 
a broader range of planetary core masses. As a measure for the 
amount of atmosphere, we use in Fig. |2 the gas pressure on the 
surface of the solid core. This avoids the ambiguity of evaluat- 
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log -^pot [erg] 


Fig. 2. Gas pressure at the bottom of the atmosphere on the surface of 
solid cores with masses between 0.1 and 5 iW® as a function of the 
nebula gas density poBC- 


Fig. 3. Density as a function of the gravitational potential for the same 
atmospheric models as in Fig.[T]for cores masses of 0.1 M® and 1 Af® 
and nebula gas densities between 5 x 10“*® g/cm^ to 5 x 10“^^ g/cm^. 


ing the mass of the atmosphere, separating the gravitationally 
accumulated gas from the background density structure of the 
circumstellar disk. Judging by the surface pressure, we see, first 
of all, that the amount of atmosphere increases with the mass of 
the planetary core. The core with 5 M® having at least 4 orders 
of magnitude more atmosphere than the 0.1 M® core. 

In Fig. 13 the findings of Fig. [1] are now placed in a broader 
scheme. Dependent on the mass of the planetary core, the re¬ 
action of the atmospheres varies between almost complete es¬ 
cape and comparatively small readjustments. For the core with 
0.1 M®, the surface pressure reduces by more than 12 orders of 
magnitude, for the 1 M® core by a factor of about 9, and by 
a factor smaller than 2 for the 5 M® core. In this simple, sta¬ 
tionary analysis, planetary cores with masses smaller than about 
0.5 M® seem unable to retain a substantial part of the nebula- 
accumulated atmosphere once the nebula disappears. Similar se¬ 
quences of atmosphere models f or a 1 M® core have been p re- 
sented by iNakazaw et al.l (Il985h and llkorna & Gendal (1200^ . 

Atmospheres around more massive cores are less sensitive to 
variations of the outer boundary density. This comparative insen¬ 
sitiveness on the outer boundary conditions is an effect which, in 
similar form, has also been observed in earlier studie s of plane¬ 
tary atmospheres with radiative oute r envelopes (e.g. JStevensonl 
[l90 IWucht^H^ ^afi^l200^ . 

The general outline of the results does not depend on the 
particular choice of physical and numerical parameters, such as 
host star and disk properties, planetary orbit, dust opacities, or 
accretion rate of planetesimals. Different parameters basically 
shift the scenario in the Mpi - Lp\ plane but leave it otherwise 
unchanged qualitatively. 

To understand the structures and properties of the atmo¬ 
spheres shown in Fig. [T] we first note that self-gravity is not im¬ 
portant in these low-mass atmospheres. The atmospheres there¬ 
fore can be considered to reside in a point-source gravitational 
potential in good approximation. In this picture, the atmospheric 
structure is independent from the particulars of the solid core 
and the surface radius of the solid core essentially becomes an 
arbitrary radius where the atmospheric structure is truncated. 

Another important observation is that, with respect to strati¬ 
fication, the atmospheres can be separated into three distinctive 
parts: ( 1 ) the inner optical thick region where convective energy 
transport takes place and the adiabatic stratification of which can 
be approximated by a polytropic structure with a polytropic in¬ 
dex of n ^ 3.3; (2) an intermediate radiative region where the 
stratification is determined by radiative transport; and (3) the 


outer part of the atmosphere that is optically thin and therefore 
accurately described by an isothermal structure (i.e., n = oo). 

The total opacity in the intermediate radiative region is usu¬ 
ally dominated by the dust opacity, which varies in a series 
of dust layers, corresponding to the formation and evaporation 
of dust species. This leads, on the one hand, to a substantial, 
nonzero temperature gradient, varying according to the local 
opacity. On the other hand, the maximum value of the tempera¬ 
ture gradient in the high-opacity regions, is limited at the adia¬ 
batic value by the onset of convective transport. On average, the 
temperature gradient in the intermediate radiative region is there¬ 
fore, somewhat smaller than, but not far off the adiabatic temper¬ 
ature gradient. As compared to the outer, optically thin part of 
the envelope, the stratification in the intermediate radiative part 
is quite similar to the inner convective part and in the tempera¬ 
ture profiles, shown in Fig. |5] and Fig. | 6 ] the boundaries of the 
convective regions are essentially indiscernible. It is therefore 
reasonable to combine both the inner convective region and the 
optical thick radiative region above in one approximated poly¬ 
tropic structure with a common polytropic index. That way, the 
planetary envelope is dissected into a two-part analytical model 
where the inner optical thick part of the atmosp here is described 
as a polytropic structure fe.g.. lEddingtonl[r926h 


p(r) = 


Cl + C2 



(9) 


while the density of the optical thin and isothermal outer part can 
be written as 

p(r) = C3expJ^^j. (10) 

Here, n is the polytropic index of the inner part, T the gas tem¬ 
perature of the isothermal outer part, and Ci, C 2 , and C 3 are 
constants. 

In this description, the density is solely a function of the 
gravitational potential. Figured which plots the same sequences 
of atmospheric structures as in Fig. [T] against the gravitational 
potential, confirms that this also holds for the numerical mod¬ 
els in good approximation. Even though important physical in¬ 
gredients of the numerical models, such as accretion luminosity, 
convection model, and realistic data for the equation of state and 
opacities, do not scale in a simple manner, the similarity between 
the two sets of solutions is very clearly apparent. 

Erom this perspective, it appears that the atmospheric struc¬ 
tures around cores of different mass are essentially homologous. 
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but shifted in the radius coordinate by the different core mass and 
trancated at different relative positions according to the radius of 
the respective core. 

The different dependency of the atmospheric structures on 
the outer boundary density can also be explained from the simple 
two-part analytical model: starting at a outer boundary with a 
density poec, a temperature Tq^c, and a radius rose, we can 
calculate the constant C 3 in Eq. [10] and obtain for the density 
profile in the outer isothermal part 


P 


= POBC exp 


1 1 

GMpi 

GMpij) 

I'^Tobc 

r 

^OBC J/ 


( 11 ) 


With increasing density toward the surface of the planet, the at¬ 
mosphere gets more opaque and at a density pxrans the isothermal 
structure changes to a polytropic stratification. In the context of 
the two-part analytical model, pxrans has to be specified as a pa¬ 
rameter. The effect of pxrans is that it scales the density in the 
deep atmosphere and thus the total mass of the atmosphere. 

The gravitational potential at this transition point 
(GMpi/rxrans) IS a function of the outer boundary density 
POBC- The temperature at the transition point Txrans is that of 
the isothermal part, i.e., Txrans = Tqbc- From these interface 
values, rxrans, Pxrans, and Txrans and integrating the pressure in a 
polytropic structure, we can calculate Ci and C 2 in Eq. j^and 
rewrite the density profile for the inner polytropic part as 


P - PXrans 



1 1 
OBC n + 1 


GMpi 


GM, 


pi 


^Trans 


( 12 ) 


The only point where the outer boundary density poBC enters this 
equation is via rxrans- For large gravitational potentials ™ » 
■ 2 ^ (i.e., r « rxrans) the density converges to 


P — PXrans 



_l _1 GMpiV 

^Tobc«+1 r j 


(13) 


independent from the density at the outer boundary. Therefore, 
the more massive a planetary core is, the less dependent is its 
deep atmosphere on the outer boundary density. This also ex¬ 
plains, as the mass of a planetary atmosphere is concentrated 
at the bottom of the atmosphere, the decreasing dependency of 
the atmosphere on the outer boundary density for more massive 
planets as illustrated in Fig. |3 

Figure|4|compares the two-part analytical model for a poly¬ 
tropic index n - 3.3 and a transition density of pxmns = 
10“^ gjerv? with the numerical results. For this comparison, we 
selected the most massive planet with 5 Me from our sequence to 
illustrate the converging density profiles in the deep atmosphere. 
The transition points between the outer isothermal and the inner 
polytropic part can be read from Fig. |4|from the intersections of 
the density prohles with the horizontal line indicating the transi¬ 
tion density. 

The very simple analytical model reproduces the numerical 
models surprisingly well. Moreover, the two constitutional pa¬ 
rameters of the analytical description are closely constrained; n 
from the atmospheric composition and pxrans scales the amount 
of atmosphere, in a similar way as the luminosity does for the 
numerical models. 


3. Dynamic outflow simulations 

3.1. Method 

For the calculation of the time-dependent atmosphere models, 
the equations of radiation hydrodynamics are solved by the 



log r [cm] 

Fig. 4. Density as a function of radius for an atmosphere around a 5 Me 
core, as computed from the detailed numerical description (TAPIR) and 
from the two-part analytical model (see text). 


TAPIR-Code using an implicit time integration scheme. The dis¬ 
cretization of the physical equations is based on finite volumes 
on a staggered mesh and advective flows are calculated using the 
second order van Leer flux limiter. 

Our system of equations consists of the equation of continu¬ 
ity 

-p + V-(up) = 0, (14) 

of 

the equation of motion 

d 47r 

(pu) + V ■ (upu) -I- VP -H pV0- /crpH -I- V ■ Q = 0 , (15) 

Of c 

the equation of internal energy 


dt 


(pe) + V ■ (upe) -H P V ■ u - AnKppiJ - 5)-i- 


+ Q : Vu -I- V ■ Pconv = 0 , (16) 


the Poisson equation 

A0 = 47rGp, (17) 


the radiation energy equation (zeroth moment of radiation inten¬ 
sity) 

■^7-hV-(u7)h-cV-H-hK: Vuh-ca'pp(7-5) = 0 , (18) 

Of 

and the radiation flux equation (first moment of radiation inten¬ 
sity) 

-^H -hV-(uH)-hcV-K-hH-Vuh-c krpH = 0 . (19) 

dt 

For the meaning of the symbols see Tab.[T] In particular, S is the 
source function of radiation with S - -T^ and Q is the viscous 

TT 

pressure tensor defined as 


Q = -qiin rpcl]: (Vu -I- (Vu)’^) - ^ V ■ u 


( 20 ) 


and scaled by the viscosity parameter q\{^. The equation of state 
enters the system of radiation hydrodynamics equations via the 
gas temperature T and the gas pressure P. Gas and dust opacities 
appear both in the form of the Rosseland-mean a-r and Planck- 
mean Kp. 
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Table 2. Standard parameters for the turbulent convection model. 


Parameter 

Adopted values 

Physical meaning 

A 

f^'ML = 2 

Mixing length, computed from 

(Ts 

0.866 

Diffusion parameter for entropy perturbation 

‘2'diss 

4 

Dissipation of turbulent energy 

^^rad 

9 

Radiative losses 


The closure of the moment description of the radiation field 
(J, H, K) requires the specification of the Eddington factor fgdd = 
K/J where we use /edd = 5 , i e-, the Eddington approximation, 
throughout. 

Convective energy transport is included in TAPIR through an 
eq uation fo r the tu rbulent convective e nergy following concepts 
of iKuhfufil jl987h along the lines of IWuchterl & Eeuchtingej 
(Il998h . Ignoring nonessential additions (turbulent pressure, con¬ 
vective viscosity, diffusion of turbulent convective energy, radi¬ 
ation pressure in the source term), the model takes the form of 


^ (piu'") + V ■ (upiu'") + p^Vad VP ■ <^'u') + 


<^diss P 


H--(21) 


A cpKr^ph? ^ 


The correlation between entropy perturbation and velocity per¬ 
turbation is provided by the diffusion ansatz 


The physical scenario of the time-dependent simulations is 
the same as for the stationary models: i.e., the computational 
domain spans from the surface of the planetary core up to the 
Hill radius where the outer boundary conditions are applied , in 
particularp = 5 x 10 g/cm^, T - 200 K, and J - S. The 
Hill radius is based on an orbit with a = lAU and a host star 
with M* = IMq. The outer boundary is open for fluid flow and 
radiation and the solution of the equation of motion and radiation 
flux equation at the outer boundary gives the mass and energy 
exchange with the disk nebula environment. 

The dust depletion factor is fixed at / = 0.01, and for the 
planetary accretion rate of planetesimals we used a constant 
mass accretion rate of Macc/ATpi = lO^^yr^^ Here, the accre¬ 
tion rate of planetesimals is considered solely as a motivation 
for the planetary luminosity, the mass of the solid core is not 
increased in consequence of accretion (and thus neither is the 
Hill radius). This inconsistency is unavoidable when simulating 
over the very long periods in time necessary to reach station¬ 
ary solutions while using the simple correlation between Macc 
and Lpi. Whereas in the stationary models, the adopted planetary 
luminosity is constant throughout the atmospheric structure, we 
assume, in the time-dependent models, that it originates at the 
bottom of the atmosphere on the surface of the solid core. 

As for the stationary models, we assumed a chemically 
homogeneous medium with sol a r com position and used gas 
opacities from lEreedman et alJ ( 2 008h. dust opa cities from 
ISemenov et al.l (12003 ). and the Saumon et al.l (Il995h equation of 
state. 


/ / ,\ A 

(i u ) ^ -PsAa/- u — , 

\ 3 or 

where the entropy gradient is evaluated as 

ds I de P dp 
dr T dr p^T dr 

and the convective flux is computed from 

Fcoav = pT <s'u') . 

Table |2] lists the standard parameters adopted for all calcula¬ 
tions, which have been adjusted so as to reproduce the so- 
lutions of the mixing lengt h theory (in the formulation of 
iKippenhahn & Weigerl]ll990l) in the stationary limit. 

Finally, the system of physica l equation can be supp lemented 
with an additional grid equation dPorfi & Drurvl[T98^ . which is 
solved together with the system of physical equations. The grid 
equation governs the motion of a fixed number of grid points 
(usually 500 points) and continuously redistributes them to al¬ 
low for an optimal resolution of the physical gradients within 
the evolving solution. 

The set of equations is solved with an implicit scheme, which 
repeatedly computes corrections to the set of variables from the 
inversion of the Jacobi-matrix until a relative accuracy of 10“^ is 
achieved. Time step control essentially aims for computational 
efficiency by balancing time step size with computational effort, 
i.e., the number of iterations required for a single step. The im¬ 
plicit time integration method allows for large time steps without 
being restricted by the Courant-Friedrichs-Lewy condition and 
therefore facilitates the efficient calculation of both fast dynami¬ 
cal processes as well as of long-term evolutionary processes. The 
simulations were terminated at an age of lO'® s as all runs are in 
a stationary state by then. To cover this period in time typically 
requires a few thousand time steps. 


( 22 ) 


(23) 

(24) 


3.2. Results 

The physical scenario of the time-dependent simulations is again 
the dynamic outflow process as described in Sect. o The evolv¬ 
ing structure of the atmospheres is shown in two examples in 
Fig.E] and Fig. | 6 ] for planetary cores with 0.1 and 2 M®, 
respectively. 

Apart from the imprints of dynamical processes, discussed 
below, the temperature profiles (lower panels) are characterized 
by an isothermal outer part and an essentially continuous tem¬ 
perature profile in the inner optical thick radiative and convective 
regions. 

In the density profiles (upper panels of Fig.|5]and Fig.| 6 ]), two 
different types of stratification are discernible: quasi-hydrostatic 
and outflow-dominated. In the former case, the stratification is 
still essentially in hydrostatic equilibrium, i.e., the density gra¬ 
dient is similar to that of the corresponding hydrostatic model. 
In the latter case, the structure is no longer in hydrostatic bal¬ 
ance but dominated by the outflow, which results in a den- 
sity gradient according to the solution of an isothermal wind 
(iLamers & Cassinefli|ll999h . In particular in the outer part of the 
upper panel of Fig. |5] the two different gradients are clearly ap¬ 
parent. 

In general, the interfaces between these two solutions would 
be discontinuous shocks, but here they are broadened by numer¬ 
ical dissipation and appear as bumps of finite width. The magni¬ 
tude of these bumps in the density profile is characterized by a 
constant total pressure (composed of gas and dynamic pressure), 
i.e.. Plot = Pgas + Pdyn - Pgas + \pu^ — coust.. The energy dissi¬ 
pation in these shocks is a source term to the internal energy and 
thus may produce spikes in the temperature profile, depending 
on the local efficiency of radiative transport. 

In Fig. |5] and Fig.| 6 ] the simulation starts from an initial hy¬ 
drostatic model (dash-dotted line) which corresponds to an at¬ 
mosphere embedded in the circumstellar disk. As soon as the 
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drop in density on the outer boundary starts, an outward flow of 
the atmospheric gas sets in. While the density profiles follow the 
general outline of the sequence of stationary models presented 
in Fig.[Tl they also show imprints of dynamical processes. First, 
during the phase of the outer boundary density reduction (i.e., the 
first 32 years) a prominent bump appears in the density profiles 
slowly traveling outward. This bump marks the interface where 
the establishing outflow from the inner region meets the quasi¬ 
hydrostatic outer part that is characterized by a smooth connec¬ 
tion to the changing outer boundary condition. In the temper¬ 
ature plot in the lower panel of Fig. |5] these shocks appear as 
spikes, in the corresponding plot for the 2 Me core in Fig. |6] 
the spikes are too small to stand out. As the outer, thin part of 
the atmosphere has a very short relaxation timescale, it can fol¬ 
low almost instantaneously the density changes imposed on the 
boundary and remains close to a hydrostatic stratification. 

In the initial expansion phase, the launch region of the out¬ 
ward flow reaches into the optically thick part where the expan¬ 
sion causes a temporary drop of the temperature below the level 
of the isothermal part, producing prominent dips in the tempera¬ 
ture profiles. These temperature dips, however, do not persist and 
are soon filled up by radiative transport. The Kelvin Helmholtz 
timescale of the relevant part of the atmospheres is about 0.3 
years for the 0.1 Me core and about 10 years for the 2 Me core, 
i.e., significantly shorter than the dynamical outflow timescale. 
At later stages in the evolution, the acceleration of the outflow 
mainly occurs in optical thin regions where radiative transport is 
much more efficient and thus prevents the formation of these low 
temperature dips. 

After the initial dynamical processes, a phase of free outflow 
follows. As the variation of the outer boundary density has al¬ 
ready reached its final (i.e., nondisk) value by then, the density 
profiles in this phase exhibit a discontinuous kink at the outer 
boundary. In this phase, the physical structure is determined 
throughout by the outflow originating in the deep atmosphere. In 
the outer part of the atmosphere, the density stratification there¬ 
fore closely resembles an isothermal wind. 

Finally, at the end of the evolution another inward travel¬ 
ing bump is produced where the dynamic outflow again meets 
an outer quasi-hydrostatic region attached to the outer boundary 
condition. This quasi-hydrostatic region slowly grows inward 
and ultimately the dynamic outflow entirely dies away and the 
whole atmosphere resettles to a hydrostatic configuration. Ac¬ 
cordingly, the final hydrostatic configuration (dashed line) again 
exactly resembles the stationary solution as illustrated in Fig. [T] 

The dynamic outflow, and in particular the bumps appearing 
between flowing and quasi-hydrostatic parts of the atmosphere, 
should not be confused with acoustic waves. Acoustic waves nat¬ 
urally appear in the solution, and the typical sound crossing time 
through the atmosphere is on the order of 10® s. Acoustic pertur¬ 
bations, however, do not contribute to mass transport and only 
very insignificantly (as compared to the luminosity) to energy 
transport. Since the resolution of acoustic waves is very much a 
numerical inconvenience in that it requires very small time steps, 
we aim to suppress acoustic perturbations in dynamic outflow 
simulations whenever possible (see Sect. l3.4l belowl. 

The timescale of the outflow and atmospheric evolution is 
studied in Fig.|7] Fig.0 and Fig.|3 The difference in timescale 
can be best appreciated from Fig.|7]showing the outflow velocity 
over the outer boundary as a function of time. For comparison, 
the speed of sound in the outer atmosphere is about 0.8 km/sec. 
The outflow velocity remains essentially constant during the out¬ 
flow phase and the duration of the dynamical outflow varies dras¬ 
tically with the mass of the solid core. Whereas the dynamical 



8.5 9.0 9.5 10.0 10.5 11.0 


log r [cm] 

Fig. 5. Density and temperature as a function of radius for the atmo¬ 
sphere around an 0.1 Me core at several snapshots during the dynamic 
outflow caused by a drop of the density at the outer boundary from 
5x10“'® g/cm^ to 5x 10“^^ g/cm^ over 32 years. The dash-dotted lines 
are the initial model, and after the dynamic outflow phase the configu¬ 
ration again settles in a hydrostatic configuration indicated by dashed 
lines. The times of the snapshot are, from top to bottom, in years: 0 
(dash-dotted line), 3.9, 7.5, 9, 11, 13, 17, 64, 141, 148, 154, oo (dashed 
line). In the color version of this figure, some lines are highlighted in 
color to simplify the identification of equivalent snapshots in the two 
panels. 


evolution last about 100 years for an 0.1 Me core, it already last 
10® years for an 0.2 M® core, 10® years for an 0.3 M® core, and 
over 10*^ years for cores more massive than 1 M®. Obviously, it 
is not sensible to interpret these long-term evolution results liter¬ 
ally as they are based on assumptions (e.g., accretion luminosity) 
that will not remain constant nearly long enough. The results are, 
however, useful to study the inherent dynamical timescale of the 
planetary atmospheres. 

The evolution of the total atmospheric mass is illustrated in 
Fig.i A more appropriate measure for the amount of atmo¬ 
sphere is, as mentioned above, the pressure at the bottom of the 
atmosphere, which is presented in Fig. |8] 

The magnitude of the evolution timescales for cores of dif¬ 
ferent mass can be verified from a simple argument based on sta¬ 
tionary models. First, we note that the density in the outer part of 
the atmosphere (not necessarily the density boundary condition) 
evolves over time in a continuous way. Therefore, the evolution 
of the atmospheric structure can reasonably be approximated by 
a sequence of stationary models with different outer densities 
such as presented in Fig. [T] The models in the sequence possess 
different atmospheric masses, and in an assumed evolution this 
mass has to be transported over the outer boundary. The density 
at the outer boundary is known from the boundary condition, 
and for the flow velocity we can put for the kinetic energy in an 
isothermal outflow Ckin = — P thus estimate u ^ Cs- 
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Fig. 6. Same as Fig. |5] except for a core with 2 A/®. The times of the 
snapshot are, from top to bottom, in years: 0 (dash-dotted line), 3.1, 3.7, 
4.4, 6.4, 2.3xl0^ 1.3xl0^ 1.6xl0«, 1.9 x lO'^, 2.9 x 10'^ oo (dashed 
line). In the color version of this hgure, some lines are highlighted in 
color to simplify the identification of equivalent snapshots in the two 
panels. 


Now provided with a difference in mass, a cross section 
(Hill sphere), a flow density, and flow velocity, we can ob¬ 
tain an estimate for the evolution time to get from one station¬ 
ary model to the other. Summing up for a sequence of mod¬ 
els yields a pseudo time-evolution of the planetary atmosphere. 
These pseudo time-evolution scales for the atmospheric mass are 
illustrated in Fig. [TO] for planetary cores between 0.1 Me and 
5 Mg). Comparison with the equivalent numerical results in Fig.|9] 
reveals that the simple flow argument is able to confirm the range 
of evolution timescales in the hydrodynamics simulations. The 
timescales derived from the sequence of stationary models are 
systematically shorter than the timescales in the dynamic mod¬ 
els. This, however, does not affect the basic argumentation but 
only reflects the simplifications of the method. While the simple 
mass flow estimate may give about the right order of magnitude, 
it cannot be expected to reproduce the results of hydrodynamical 
simulations with any accuracy. 

The stark contrast in the atmospheric reaction (measured in 
Fsurf) to the drop in outer boundary density between cores of 
different mass has already been illustrated in Fig |2] In addi¬ 
tion, Fig.[8]also highlights the very different timescales on which 
these reactions happen: Therefore, not only are the consequences 
for the atmospheric structure to the removal of the embedding 
nebula gas for a massive core much smaller than for a low mass 
core, the dynamic reaction is also much slower and thus the at¬ 
mosphere of the massive core is unlikely to evolve very far in 
any reasonable time. 



Fig. 7. Outflow velocity over the outer boundary as a function of time 
in dynamic outflow simulations for atmospheres around planetary cores 
of different mass. 



log t [yr] 


Fig. 8. Atmospheric pressure at the surface of the solid core as a func¬ 
tion of time in the dynamic outflow simulations. The line-type coding 
is the same as in Fig.[7l 



Fig. 9. Total mass of the atmosphere in the dynamic outflow calcula¬ 
tions as a function of time for planetary cores of various masses. The 
line-type coding is the same as in Fig.[7l The different initial curvature 
of the lines, as compared to Fig. [8] is due to the fact that the total mass 
not only includes gravitationally accumulated gas but also background 
disk structure within the Hill sphere that can contribute significantly in 
the case of thin atmospheres. 


3.3. Physical parameters 

To obtain a measure for the dependencies of the results on both 
the physical and numerical parameters involved in the calcula¬ 
tions, we investigate here, and in Sect. 13.41 a parameter space in 
the most obvious and important parameters. 
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Fig. 10. Total atmospheric mass of stationary models around cores 
0.1 and 5 plotted on a time axis calculated from a simple model 
for the transition between individual stationary atmospheres (see text). 


Apart from the mass of the solid core, the two most impor¬ 
tant constitutive parameters of stationary models of planetary at¬ 
mospheres are the planetary luminosity, here associated with a 
mass accretion rate, and the dust grain depletion factor /. Un¬ 
fortunately, both of these parameters are very poorly constrained 
both from theory and observations. The effect of these param¬ 
eters is illustrated in the first and second panel of Fig. [TT] and 
Fig. [12] for solid cores of 0.1 and 2 Me, respectively. 

The most obvious effect of both the planetary luminosity and 
the dust depletion factor is the scaling of the amount of atmo¬ 
sphere in the initial models and thus of the atmospheric surface 
pressure. The sensitivity to the parameters is about the same in 
the two examples, even though the difference is harder to pick 
out in case of 0.1 Me core because of the larger scale of the 
pressure axis. 

The influence of the outer boundary conditions on the atmo¬ 
spheric models has been tested by varying the orbital distance 
and the magnitude of the drop in nebula density in the third and 
fourth panel of Fig. [TTjand Fig. [ 12 ] 

The orbit of the planet enters Eq. [T]for Hill radius and thus 
defines the localization of the outer boundary of the model. As 
both the orbital radius and the mass of the host star only enter in 
this equation, it is sufficient to consider here the orbital distance 
only. We do not take the effect that a different orbit or host star 
mass will also imply a different disk gas temperature into ac¬ 
count. Figure [TT] and Fig. [ 12 ] show, however, that the placement 
of the outer boundary, as specified through the orbital distance, 
does not appreciably affect the results. This is reassuring insofar 
as the placement of the outer boundary is a somewhat arbitrary 
decision. 

The magnitude of the density drop, i.e., the assumed remain¬ 
ing ambient density after the protoplanetary disk has evaporated, 
only has the effect of defining the overall duration of the outflow 
phase. The physical behavior and the intrinsic timescale are un¬ 
changed. Thus, the simulations evolve along a very similar pat¬ 
tern but are terminated at a stage according to the final ambient 
density value. In terms of outflow duration, the effect is larger for 
the 2 M® core than for the 0.1 M® core because of the different 
slope of the atmospheric evolution tracks. 

Finally, in Fig. [13] we illustrate the effect of the tempera¬ 
ture of the disk environment on the atmospheric structure and 
evolution. In general, the lower the temperature of the embed¬ 
ding nebula, the more dense (and thus the more massive) is a 
stationary envelope. For the two-part analytical model, this is 
directly evident from Eq. [TT] and Eq. [12] As the amount of at- 



Fig. 11. Effect of important physical parameters on the surface gas pres¬ 
sure as a function of time in dynamical simulations of an atmosphere 
around an 0.1 Me core. 

Erom top to bottom: Accretion rate of planetesimals onto the 
planetary core (defining the planetary luminosity); dust grain 
depletion factor /; orbital radius of the planet; magnitude of 
the drop in nebula gas density, starting in each case from 5 x 
10 “'° g/cm^. 

mosphere is a main factor determining the time evolution in the 
outflow scenario, the model series with different temperatures in 
Eig.[l3]looks strikingly similar to a series with varying core mass 
(Fig .[ 8 ]l. The large spread in evolution timescales in Eig.[T3] how¬ 
ever, is much reduced when focusing on planets in the habitable 
zone as this imposes strict limits on the ambient temperature. 

3.4. Numerical parameters 

Apart from physical parameters discussed above in Sect. 13.31 
several purely numerical parameters affect the solution and de¬ 
serve discussion. 

Convection parameters: The turbulent convection model in¬ 
volves several parameters, most importantly the mixing length 
parameter oml- The set of standard parameters, adjusted with 
reference to the mixing length theory of convection, has been 
compiled in Tab. [2] Overall, the effect of the convection parame¬ 
ters remains small. Eor reasonable choices of auL, the structure 
is very close to the adiabatic one at all times. Therefore, the tradi¬ 
tional and more simple approach to the modeling of convection 
zones in planetary atmospheres of directly using the adiabatic 
gradient once the Schwarzschild criterion is fulfilled is well jus- 
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Fig. 12. Same as Fig. [n for a planetary core of 2 M®. 
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Fig. 13. Atmospheric pressure at the surface of a 0.3 core as a 
function of time, calculated for different temperatures of the surround¬ 
ing nebula. 

tified. In the present study, however, this sort of approach would 
not be applicable as hydrodynamic simulations also require a 
dynamical convection model. 

Viscosity: Dissipation enters the scheme in several distinctive 
ways: from temporal discretization, spatial discretization, and 
from artificial numerical viscosity. One particular problem in the 
dynamical escape simulations is the need to suppress acoustic 
waves, which would travel back and forth in the thin outer atmo¬ 
sphere as soon as the density variation at the outer boundary be¬ 
gins. These acoustic waves do not contribute to the evolution of 
the atmospheric structure, and yet the small time steps necessary 


to resolve these waves would make it impossible to cover the 
relevant evolutionary timescales in all but a few cases. Acoustic 
perturbations in planetary atmospheres are subject to physical 
dissipation processes, such as radiative cooling, turbulence, and 
non-radial dispersion. These effects, however, are either small 
(radiative cooling) or not included in the model. Therefore, de¬ 
pending on the particularities of the atmospheric model, some¬ 
times significant numerical viscosities had to be deployed to 
raise the overall dissipation of the code to levels high enough 
to suppress such acoustic waves. 

So while some dissipation remains indispensable, it is pos¬ 
sible to disentangle the effect of the individual contributions to 
some degree. 

First, temporal discretization errors do not affect the solu¬ 
tions spatially but only dissipate fast (high frequency) physical 
variations. In the application to evolution simulations, this is a 
very desirable effect and thus we used a first order time inte¬ 
gration scheme for all time-dependent calculations. Dissipation 
from temporal discretization errors is the most important means 
to suppress acoustic perturbations and therefore a second order 
scheme would only be advisable when one is interested in oscil¬ 
latory phenomena. 

The effect of the artificial viscosity can be studied from the 
comparison of simulation runs without (for a model where this is 
possible) and with various amounts of viscosity. In the produc¬ 
tion runs viscosity parameters between 0 and 1 have been used, 
and in this range the viscosity only moderately affects the over¬ 
all results, the largest effect being apparent for the least massive 
cores. But even in the case of the 0.1 M® core, the effect is only 
a slight change in the exponential slopes in the evolution of vari¬ 
ables with a corresponding change in evolution time, the solution 
otherwise remains very similar qualitatively. 

Dissipation from spatial discretization scales with the spatial 
resolution to the same order as the accuracy of the numerical 
scheme, which here is essentially of second order. To estimate 
this effect, we changed the number of grid points from the de¬ 
fault number of 500 in a range between 100 and 5000. Another 
way to reduce spatial discretization errors is to enable the adap¬ 
tive grid equation, which continuously redistributes a fixed num¬ 
ber of grid points according to gradients in the evolving solution 
and thus increases the effective spatial resolution of the physical 
variables. 

The result of both these numerical tests is that the effect 
of the spatial discretization errors on the solutions of the time- 
dependent calculations is comparable to, but smaller than, the ef¬ 
fect of artificial viscosity (in the adopted parameter range). From 
this we can estimate that, for the results presented in Sect. 13.21 
the inaccuracies introduced from the limited spatial resolution 
are smaller than the (controlled) effect of the artificial numerical 
viscosity. 

4. Discussion 

For the astrophysical interpretation of the numerical results, it is 
important to keep in mind the assumptions and limitations in¬ 
herent to the numerical scheme. First, the assumption of a con¬ 
stant planetary luminosity during the whole evolution scenario is 
clearly not realistic. In the disk-phase, ga s-drag is a possible way 
to mediate migration of planetesimals dWeidenschillind fl 97711 
and thus to provide a resupply for the accretion of planetesimals 
onto the planet. It seems likely that this changes fundamentally 
once the disk evaporates and in later stages of the evolution of the 
planetary system, the accretion rate will soon become insignifi¬ 
cant for the energy balance of the planetary core. 
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Second, the numerical treatment of some atmospheric 
physics is rather simple; chemical inhomogeneities are neglected 
and the simple dust model negates effects such as grain growth or 
atmospheric rain-out. The radiation field, which, cotTesponding 
to the complex materials physics, will have a very feature-rich 
spectral distribution, is approximated by gray radiative transport 
and frequency averaged opacities. Furthermore radiative trans¬ 
port is essentially diffusion-like and iri'adiation from the host star 
onto the planetary atmosphere is not explicitly modeled. 

While some of the numerical limitations might be lifted by 
future advances of the code, some other restrictions, e.g., the un¬ 
certainties in the modeling of dust, seem unavoidable. 

However, it is not the scope of this paper to predict realistic 
evolution scenarios of specific planets, nor to construct detailed 
structure models for planetary atmospheres, but to explore the 
general trends and the dynamical behavior of disk-accumulated 
planetary protoatmospheres. As the test-calculations against 
physical parameters presented in Fig. [TT] and Fig. [12] illustrate, 
the general physical behavior is not qualitatively affected, even 
if one allows for significant variance in important parameters. In 
particular, there appear to be some robust conclusions that can 
be drawn from the numerical results. 

The stationary atmosphere models demonstrate, that the de¬ 
pendency of the atmospheric structure on the sutTounding en¬ 
vironment varies in a very nonlinear way with the mass of the 
solid core. There is a distinct cutoff between low-mass cores, 
which lose almost all of their atmosphere once they are removed 
from the disk nebula environment and more massive cores that 
keep a substantial atmosphere. The position of this cutoff in core 
mass, however, is subject to physical parameters that both reflect 
uncertainties in modeling as well as actual differences from one 
planetary system to the other. The analytic approximation re¬ 
veals that this difference corresponds to atmospheres with deep 
adiabatic parts that are comparatively independent of the outer 
boundary conditions and atmospheres around small cores, which 
are almost entirely isothermal and strongly scale with the density 
at the outer boundary. 

The time-dependent simulations exhibit an equally drastic 
difference in the timescale of atmospheric outflow. Low-mass 
planetary cores not only lose almost all of their atmosphere 
with the evaporation of the circumstellar disk, they also do so 
very quickly compared to other evolutionary processes. More 
massive cores, on the other hand, show very long evolution 
timescales and it seems clear that atmospheric mass loss by 
dynamic outflow is not efficient for massive cores, and con¬ 
sequently they will retain a large part of the atmosphere ac¬ 
creted in the disk-embedded phase. The cooling of the planet in 
later stages of evolution and the decreasing planetary luminos¬ 
ity will only enhance this trend. Hydrodynamic atmospheric loss 
caused by the high XUV radiation of the protoplanet’s young 
host star is also unlikely to remove a significant fraction of 
nebula-capt ured hydrogen envelopes around more ma ssive plan¬ 
etary cores (lErkaev et alj|201^ iLammer et al.ll2014l) in orbital 
distances related to the habitable zone of solar-like stars. More¬ 
over, nonthermal atmospheric escape processes, such as pick up 
of planetary hydrogen ions by the stellar wind, are also too weak 
to remove dense nebula-captured hydrogen envelopes from so- 
called “super-Earths” w i th cor e masses above 1.5 Earth-masses 
(iKislvakovaetalJ 12013112014 at 1 AU. As we have shown in 
Sect. 13.21 the range of evolution timescales in the hydrodynam¬ 
ics simulations can also be understood from a simple model for 
the transition between individual models in a sequence of sta¬ 
tionary atmospheres. 


A very general conclusion one can draw from the time- 
dependent simulations is that the atmospheric relaxation 
timescales frequently become very long and that one cannot rely 
on the assumption that a planetary atmosphere will always be 
close to the stationary solution. In particular, the atmospheric 
relaxation timescales can also become very large with respect 
to the evolution and evaporation timescale of the protoplanetary 
disk. Therefore, numerical modeling of disk-accreted planetary 
atmospheres should not only rely on stationary models but be 
based on time-dependent methods. 
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